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ABSTRACT 

We report on the results of the first 3D SPH simulation of massive, gravitationally unstable proto- 
planetary disks with radiative transfer. We adopt a flux-limited diffusion scheme justified by the high 
opacity of most of the disk. The optically thin surface of the disk cools as a blackbody. The disks 
grow slowly in mass starting from a Toomre-stable initial condition to the point at which they become 
marginally unstable. We find that gravitationally bound clumps with masses close to a Jupiter mass can 
arise. Fragmentation appears to be driven by vertical convective-like motions capable of transporting 
the heat from the disk midplanc to its surface on a timescale of only about 40 years at 10 AU. A larger or 
smaller cooling efficiency of the disk at the optically thin surface can promote or stifle fragmentation by 
affecting the vertical temperature profile, which determines whether convection can happen or not, and 
by regulating the accretion flow from optically thin regions towards overdense regions. We also find that 
the chances of fragmentation increase for a higher mean molecular weight, (i, since compressional heating 
is reduced. Around a star with mass 1-M© only disks with /i > 2.4, as expected for gas with a metallicity 
comparable to solar or higher, fragment. This suggests that disk instability, like core-accretion, should be 
more effective in forming gas giants at higher gas metallicities, consistent with the observed correlation 
between metallicity of the planet-hosting stars and frequency of planets. 

Subject headings: accretion disks - - hydrodynamics - - planetary systems:formation — solar sy 
stem:formation — methods:N-Body simulations 



1. INTRODUCTION 

Long-lived clumps can form in a massive gravitational 
unstable protoplanetary disk if the gas is isothermal 
(Mayer et al. 2002) or cools on a timescale comparable 
to the local orbital time (Rice et al. 2003a, b, Mayer et al. 
2004a, Rice et al. 2005, Mayer et al. 2005, Mejia et al. 
2005). Fragmentation under such favourable conditions 
is a numerically robust result obtained with SPH codes 
as well as grid-based codes with both static and adaptive 
meshes (Durisen et al. 2005, Mayer et al., in prep.). The 
main open question is whether the required short cooling 
times can be obtained in the disks. Shock heating along 
spiral arms is intense during the phase of strong gravita- 
tional instability (Pickett et al. 2000, 2003; Mayer et al. 
2004b) and will tend to erase clumps exactly when they 
have the best chance to collapse. Boss (2002a,2002b,2004) 
has used a grid-based code which solves radiation trans- 
port with the diffusion approximation. He found that the 
disk cools rapidly through convection instead of radiation 
and forms long lasting clumps. Cai et al. (2006) adopt 
flux-limited diffusion in their cylindrical grid code. They 
find long cooling times, no evidence of convection and no 
fragmentation. Their method is very different from that 
of Boss (2004). While Boss does not use a flux-limiter but 
embeds the disk in a thermal bath at a constant tempera- 
ture in the range 40-50 K, Cai et al add an Eddington-like 
atmosphere on top of the optically thick layer, explicitly 
matching the fluxes at the boundary. Nelson et al. (2000) 
also did not find fragmentation in a 2D SPH simulation 



to which a plane-parallel atmosphere radiating as a black- 
body was added assuming that rapid vertical cooling was 
achieved through convection. Here we present the first 3D 
SPH simulations of protoplanetary disks with flux-limited 
diffusion. 

2. THE SIMULATIONS 

The disk models are set up as described in Mayer et 
al. (2004). Initially disks extend from 4 to 20 AU, have 
a power-law surface density profile and have a minimum 
temperature of 40 K at the outermost radius. The simula- 
tions employ 10 6 gas particles with a gravitational soften- 
ing 0.06 AU, while the central star is represent by a particle 
with mass 1M Q and softening 0.2 AU. Disks are grown in 
mass from a stable state with a high Toomre parameter 
(Tommre 1964), Q m i n > 4. Mass growth is stopped ei- 
ther when fragmentation begins or when the disk reaches 
a mass equal to 0.2M Q . Thanks to the very high resolution 
we always fulfill the criterion of Bate & Burkert (1997) to 
avoid spurious fragmentation. The simulations were run 
with the parallel SPH+N-Body code GASOLINE (Wads- 
ley, Stadel & Quinn 2004). We solve the energy equation in 
asymmetric form accounting for irreversible shock heating 
via the standard Monaghan artificial viscosity term (we 
include the Balsara correction term to reduce unwanted 
shear viscosity). 

The radiation transport is implemented using the dif- 
fusion approximation and the flux-limiter of Bodenheimer 
et al. (1990). A similar method has been used by White- 
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house & Bate (2006) to study the collapse of molecular 
cloud cores. The flux limiter appears as a coefficient of a 
diffusive term in the energy equation. Following Cleary & 
Monaghan (1999), the energy equation reads 



Am b k a k b 
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where the summation is over neighboring particles, W is 
the smoothing kernel, r a b the vector from the position of 
particle a to particle b, and k a takes the form of a thermal 
conductivity term 
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where n a is the opacity, a is the Stefan-Boltzmann con- 
stant and A a is the flux limiter. This is a stable method 
that relaxes the noise of second derivatives. 




FlG. 1. — Color coded logarithmic density plots of the disk 
seen face-on after 1300 years of evolution. From top left to 
bottom right we show the runs with (a) /i = 2.4, EDA (edge 
detection angle) =30 degrees, (b) /j, — 3, EDA=60 degrees, (c) 
fi = 2.4, EDA=40 degrees, (d)fi = 2.7, EDA=40 degrees. The 
maximum densities shown, at the clumps position, reach values 
higher than 10 _9 g/cm 3 . 

To implement the above expression in Gasoline, we 
add the calculation of the temperature gradient to the 
smoothing operation already in place to calculate the force 
due to the pressure gradient. The (Rosseland mean) opac- 
ity of each particle is interpolated from a table of n(T,p) 
(D'Alessio et al. 1997). 

The next step is to allow particles on the boundary of 
the disk to radiate their energy away to infinity. To find 
particles that are "on the edge" of the disk, we examine the 
directions to all of the neighbors used in smoothing sums. 
If a particle has no neighbors within a certain fraction of a 
solid angle from a preferred direction (edge detection an- 
gle, hereafter EDA) it is considered an edge particle. From 
the geometry of a disk, the preferred directions (treated in- 
dependently) are out of the plane of the disk (both up and 



down) and radially outward. We let edge particles radiate 
as blackbodies, adding a term to the energy equation of 
each particle, 

U a = f a SaT*/m a . (3) 

where S — iirh^ is the surface through which the particle 
radiates and h a is the smoothing length of the particle. 
The "edgeness factor" f a represents the fraction of their 
surface area over which a particle radiates. It is usually 
zero, and takes value 1/2 for particles on one of the up, 
down, or out boundaries, 1 for those on the edge in both 
the up and down directions, and 3/4 for those on the edge 
in the out and either up or down direction. 

The EDA used in the check for "edgeness" is formally a 
free parameter. We only consider angles in the range 30-60 
degrees (measured as linear angle in a cone centered on a 
given particle) since with these the particles identified as 
"edge" do lie in regions where the optical depth is r < 1. 
Different choices for the EDA will effectively yield a dif- 
ferent radiative efficiency since both the radiating surface 
area and the edgeness factor for each particle will change. 
The smallest possible size of the radiating surface area is 
that of the geometric surface area of the disk, which cor- 
responds to an EDA of 50 degrees when the disk enters 
the phase of gravitational instability. Within the range 
of angles considered the radiative surface area changes by 
about a factor of 2. 

The term in equation 3 has a timescale that can be sev- 
eral orders of magnitude smaller than the others. To han- 
dle this, we apply a special Burlisch-Stocr stiff integrator 
to this term over a single timestep in the outer leapfrog in- 
tegrator. We hold the other terms of the energy equation 
constant over this short integration. 

A near solar opacity is used in most of the simulations. 
The (mean) molecular weight of the gas is by default equal 
to the solar metallicity value (p = 2.4) but is varied be- 
tween /j = 2 (pure molecular hydrogen) and \x = 3 (highly 
metal rich). We have run about 15 simulations with vary- 
ing molecular weight and opacity. 

3. COOLING, HEATING AND FRAGMENTATION 

The disk grows uniformly in mass at the rate of ~ 
lO _4 M /yr, approaching 0.1M Q after about 10 3 years. 
When its mass grows above 0.05M Q the Toomre param- 
eter Q drops below 2 in the outer part of the disk and 
strong spiral patterns begin to appear. The shocks occur- 
ring along the spiral arms limit the growth of their am- 
plitude as the increasing pressure counteracts self-gravity. 
Yet, fragmentation occurs in some of the simulations once 
the disk mass is in the range 0.12 — 0.15Af Q (Figure 1). 
Whether the disk fragments or not depends on the details 
of thermodynamics in the disks, in particular the molecu- 
lar weight of the gas and the extent of the radiative layer 
(see Figure 1). 

We find that a molecular weight equal to 2.4 or larger 
is in general necessary for fragmentation to occur (Figure 
1). Larger molecular weights effectively soften the pressure 
gradients across the spiral arms and allow fragmentation 
to happen for a smaller radiating surface area (Figure 1). 
This sensitivity is reminiscent of the sensitivity on the adi- 
abatic index 7 found by Mayer et al. (2004b) and Rice et 
al. (2005) (here we fix 7 = 7/5). Conversely we varied 
the disk opacity by a factor of 50 and found almost no 
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difference in the outcome, confirming the results of Boss 
(2002a). Indeed when cooling at the surface is switched 
off the temperature in the midplane is the same as in a 
control run performed with an adiabatic equation of state. 
This is not surprising. The timescale over which thermal 
energy is transported from the midplane to the surface, 
which is of order 10 4 years (Boss 2003), is the product 
of the radiative diffusion time, which depends on opacity, 
times the ratio between gas pressure and radiation pres- 
sure. The latter factor is orders of magnitude higher than 
the radiative diffusion time, is independent of opacity, and 
is what determines the radiative cooling time of most of 
the disk. This confirms earlier findings by Pickett ct al. 
(2000, 2003) on the crucial role of compressional/shock 
heating in disk instability. 

However, the cooling time has to be ~ 100 times shorter, 
namely of order of the orbital time, for fragmentation to 
happen (Rice et al. 2003a; Mayer et al. 2004a). Such 
short radiative cooling times are only possible at the op- 
tically thin edge in our disks. Boss (2002b, 2004) finds 
that the disk can cool on such short timescales thanks to 
convection. The optical depths in our disks are very large, 
t > 10, near the disk midplane. Convection is expected 
(Ruden & Pollack 1991) in such conditions. The follow- 
ing sequence of events is observed in the simulations and 
translates into the temperature evolution of overdensities 
shown in Figure 2; first, density maxima develop along an 
optically thick spiral arm in the disk midplane; second, 
the gas in the spiral arm is shock heated to temperatures 
above 150 K; third, if the disk atmosphere is cold enough 
while the midplane is shock heated vertical upwcllings and 
downwellings begin which cool the clumps and allow them 
to become gravitationally bound. Such vertical motions 
have typical speeds of order 0.1-0.4 km/s (the typical or- 
bital speed at the same radii is ~ 10 km/s). They involve 
most of the gas particles in a locally overdense region of 
the disk and occur when the vertical entropy gradient is 
superadiabatic. These motions where not observed in pre- 
vious calculations with fixed equations of state in which 
steep vertical temperature gradients were not produced 
(Mayer et al. 2004b). The pressure scale height, whose 
scale should be comparable to the mixing length in the 
convective motions, is from 10 to 50 times bigger than the 
SPH smoothing length, hence convection should be well 
resolved. 

The measured vertical speeds are high enough to re- 
distribute thermal energy on a timescale of about 30-50 
years, i.e. from comparable to slightly longer than the 
orbital time. With such short cooling timescales fragmen- 
tation is expected (Mayer et al. 2004a; Rice et al. 2003a). 
It is difficult to measure the net transport of energy due 
to convection because turbulent motions can arise also as 
a result of shocks (Cai et al. 2006) and simultaneously 
vertical motions are also produced by the accretion flow 
towards the overdensities at the midplane. 
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FlG. 2. — Top: Evolution of the vertical temperature pro- 
file of an overdense region that collapses and becomes a grav- 
itationally bound clump in the run with /i = 2.7 and EDA 
= 40 degrees (same clump as in Figure 3). The profiles are 
taken at t= 1220 years (solid line), t= 1221years (dashed line), 
t= 1222 years (short dot-dashed line) and t= 1223 years (long 
dot-dashed line), and are azimuthally averaged. As the pro- 
file becomes steeper the gas becomes convectively unstable. It 
then relaxes to a shallower, more stable profile, with an av- 
erage temperature lower than at the start (dot-dashed line). 
Bottom: Velocity component perpendicular to the disk mid- 
plane as a function of the distance from the center of mass of 
the clump at t= 1221 years. 

In addition to affecting convection by changing the ver- 
tical temperature profile, the radiative efficiency at the 
cold optically thin layer has a second influence on the final 
outcome of the disk instability. The dense regions along 
the spiral shocks can engulf material from colder regions 
outside the shock front (Figure 3), and this accretion flow 
is stronger for colder material. These regions are optically 
thin down to a vertical scale height half of that of the 
neighboring overdense spiral arms and are naturally pro- 
duced as the arms unwind. They provide a reservoir of cold 
gas that overdensities can accrete, promoting their growth. 
This is an example of the complex coupling between the 
dynamics and thermodynamics seen in these simulations. 
Nevertheless, we notice that the typical midplane temper- 
atures (> 200K) and midplane surface densities (~ 2000 
g/cm 2 ) of our disks where clumps first appear (at about 
15 AU) are admitted as fragmenting solutions by the ana- 
lytical model of Rafikov (2005), given a Toomre parameter 
Q = 1.4 and a cooling time t coo i — 1.5T or & as the criteria 
for fragmentation (see Mayer et al. (2004a). 

Once formed, gravitationally bound clumps have masses 
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that range range from one to a few Jupiter masses, are 
differentially rotating and have densities a million times 
higher than the background density. They reach r ~ 100 
at the scale of the softening and temperatures T > 300 
K. In general they resemble clumps formed in simulations 
with an adiabatic equation of state above a fixed density 
threshold (Mayer et al. 2004a,b). 

4. DISCUSSION 

We have shown that disks can fragment even when ra- 
diative transfer is included. However, fragmentation re- 
quires a disk somewhat more massive than in previous 
works relying on simple equations of states or parameter- 
izations of the cooling time., M > 0.12M Q as opposed to 
~ Q.IMq (Mayer et al. 2004a,b). It is strongly dependent 
on the efficiency of cooling at the optically thin boundary 
and on the molecular weight. The newly discovered depen- 
dence on molecular weight is intriguing since it establishes 
for the first time a direct link between the disk instabil- 
ity mechanism and the metallicity of the gas. Such a link 
is straightforward in the competing model, core-accretion 
(Pollack et al. 1996). Hence the formation of giant planets 
by disk instability should be favoured around more metal 
rich stars, in line with the observed correlation (Fischer 
& Valenti 2005). A more metal rich disk will also have 
a higher cooling rate in the optically thin region (Cai et 
al. 2006), which also appears to promote gravitational in- 
stability, whilst the associated increase of opacity in the 
optically thick regions will not affect the outcome based 
on our results. 




FlG. 3. — Color coded logarithmic temperature plots of a re- 
gion of size about 4.5 AU centered on one of the gravitationally 
bound clumps in the run with \x = 2.7 and EDA=40 degrees. 
Left: a smoothed temperature map is shown (the temperature 
range goes from 10 K (dark blue) to 250 K (dark red) at time 
t = 1248 years showing the cold, low density regions adjacent 



to the hot gas in the spiral shock and clump. Right: tem- 
perature map of the gas particles surrounding the same clump 
about half an orbit earlier with the particles that will accrete 
by t = 1248 years marked in green. The brighter the color the 
higher the temperature. Most of the accreted mass appears to 
come from the colder regions outside the spiral shock. 

The radiation physics in the simulations presented here 
is still simplified. In particular, the disk edge could be 
heated as a result of irradiation from the central star 
and/or neighboring bright stars, erasing the steep temper- 
ature gradient necessary to drive convection. On the other 
end, when an overdensity forms dust might accumulate at 
its location owing to the action of gas drag (Haghighipour 
& Boss,. 2003, Durisen et al. 2005), thus increasing lo- 
cally the molecular weight and enhancing the chances of 
fragmentation. In the latter case the clumps would be 
automatically enriched in heavy elements relative to the 
surrounding disk, as seen in the gas giants of the Solar 
System (Saumon & Guillot 2004) and predicted by recent 
models of core accretion (Alibert et al. 2005). 

There are situations in which fragmentation may be 
achieved more easily than suggested here. Indeed a new set 
of simulations being completed as we write indicates that if 
the mass of the central star is lowered to M = 0.5M Q , the 
disk fragments for masses < 0.1 M for a molecular weight 
as small as [i — 2. In this case a lower shear allows a higher 
pressure front to be tolerated from growing overdensities. 
Therefore with disk instability formation of giant planets 
around M dwarfs is likely (see also Boss 2006). We also 
find that if the disk does not grow slowly in mass but 
starts with Q m in — 1-4, as it might happen in response 
to a sudden external perturbation, fragmentation occurs 
with [i — 2 and edge detection angles of up to 60 degrees. 
On the other end, previous work shows that long lasting, 
time-dependent perturbations such as those induced by 
the presence of a binary companion might quench the frag- 
mentation by increasing the heating due to spiral shocks 
(Mayer et al. 2005). 

In summary, it seems that an ultimate answer requires 
that all the details of the radiation physics in the disks and 
in their surrounding environment are taken into account 
together with the details of the disk formation process. 
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